%%%Main File for BCDR Simulations - gains from industrial policy and trade
%%%policy

%This file computes the gains from a given trade and industrial policy
%using the SOE assumption

%CES Version: ROBUSTNESS CHECKS for non-manufacturing SE and upper tier EoS
%(Figure 5 in BCDR)

clear all;

%digits(64);

K=32; %number of sectors
N=62; %includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows

%Note: data is sorted by sector, then origin.

A=importdata('simulation_ge_data.csv');
B=importdata('sim_params_no_int.csv'); %these are the estimated elasticities

%Data vectors
X_ink = A(:,6:5+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
V_ik=A(:,3); %vector of sectoral value added, (NxK)x1
expenditures_ik=A(:,5); %vector of expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1

%Parameter vectors
theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
h_theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,2);zeros(K-2-K_man,1)]; %trade elasticities
h_gamma = [zeros(2,1);B(:,1);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
gammaXtheta=h_theta.*h_gamma;
rho=0.87; %Upper tier elasticity of substitution


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%Other data matrices
V_ik = reshape(V_ik, [N 1 K]);
V_i=sum(V_ik,3);
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

%Create trade shares, expenditure shares
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

expenditures_n=sum(expenditures_nk,3);
beta_nk=expenditures_nk./repmat(expenditures_n,[1 1 K]);

%Create sectoral domestic trade shares
dom_lambda_ik=zeros(N,K);
for i=1:N
    for k=1:k
        dom_lambda_ik(i,k) = lambda_ink(i,i,k);
    end
end
        
%Create total sectoral exports
exports_ik=zeros(N,K);
for i=1:N
    for k=1:k
        exports_ik(i,k) = V_ik(i,k)- X_ink(i,i,k);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy: Different upper tier elasticities of
%substitution

%Do a range of rho, store averge gains from industrial policy only
rho_range=.37:.04:2.07;
gfoip_rho_robust=zeros(size(rho_range,2),1);
gfotp_rho_robust=zeros(size(rho_range,2),1);
gfop_hte=zeros(N-1,1);
gftp_hte=zeros(N-1,1);
gfip_hte=zeros(N-1,1);

theta=h_theta; 
gamma = h_gamma;
counter=0;

for j=1:size(rho_range,2)
    rho=rho_range(1,j);
    for i=1:N-1
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Optimal trade and industrial policy
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfop_hte(i,1)=model_output.welfare;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Trade policy only
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gftp_hte(i,1)=model_output.welfare;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Industrial policy only
    %Specifying taxes and subsidies
    t_k=zeros(K,1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfip_hte(i,1) = model_output.welfare;
    end
    
    
    %Compute average gains from industrial policy and trade policy
    gfoip_rho_robust(j,1)=sum(gfop_hte-gftp_hte)/(N-1);
    gfotp_rho_robust(j,1)=sum(gfop_hte-gfip_hte)/(N-1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy: Different non-manufacturing scale
%elasticities

%Do a range of range of non-manufacturing gamma, store averge gains from industrial policy only
gamma_nm_range=0:.005:.21;
gfoip_gamma_robust=zeros(size(gamma_nm_range,2),1);
gfotp_gamma_robust=zeros(size(gamma_nm_range,2),1);

theta=h_theta; 
rho=.87;
counter=0;

for j=1:size(gamma_nm_range,2)
    gamma=[gamma_nm_range(1,j)*ones(2,1);B(:,1);gamma_nm_range(1,j)*ones(K-2-K_man,1)];
    for i=1:N-1
    %Defining data vectors and parameters
    V = V_i(i); %total value added
    X=exports_ik(i,:)'/V; %sectoral exports (divided by value added)
    D= D_i(i)/V; %Aggregate trade balance (divided by value added)
    beta = squeeze(beta_nk(1,i,:)); %sectoral expenditure shares
    lambda= dom_lambda_ik(i,:)';  %own expenditure shares by sector
    V_k = V_ik(i,:)'/V; %value added by sector (divided by value added)
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Optimal trade and industrial policy
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfop_hte(i,1)=model_output.welfare;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Trade policy only
    %Specifying taxes and subsidies
    t_k=1./(theta+1);
    s_k=zeros(K,1);
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gftp_hte(i,1)=model_output.welfare;
    
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Industrial policy only
    %Specifying taxes and subsidies
    t_k=zeros(K,1);
    s_k=gamma;
    %Solving the model
    model_output=gfop_soe_ces_solver(V_k,D,X,lambda,beta,gamma,theta,rho,t_k,s_k,K);
    gfip_hte(i,1) = model_output.welfare;
    
    end
    
    %Compute average gains from industrial policy and trade policy
    gfoip_gamma_robust(j,1)=sum(gfop_hte-gftp_hte)/(N-1);
    gfotp_gamma_robust(j,1)=sum(gfop_hte-gfip_hte)/(N-1);
end

gains_ces_soe_rob_rho=[rho_range' gfoip_rho_robust gfotp_rho_robust];
gains_ces_soe_rob_gamma=[gamma_nm_range' gfoip_gamma_robust gfotp_gamma_robust];

csvwrite('gains_ces_soe_rob_rho.csv',gains_ces_soe_rob_rho);
csvwrite('gains_ces_soe_rob_gamma.csv',gains_ces_soe_rob_gamma);
